前言

………

这篇博客作为(introduction to optimization)这篇笔记文的后续

线搜索

线搜索首先需要进行划界(bracket)

有没有什么办法可以加快划界的速度呢？有，逆抛物内插（Inverse Parabolic Interpolation）就是一种技术，它可以使得划界算法超线性收敛(brent方法)。

总的来说, 一个线性搜索过程如下

一般使用brent方法

In numerical analysis, Brent’s method is a root-finding algorithm combining the bisection method, the secant method and inverse quadratic interpolation. It has the reliability of bisection but it can be as quick as some of the less-reliable methods. The algorithm tries to use the potentially fast-converging secant method or inverse quadratic interpolation if possible, but it falls back to the more robust bisection method if necessary. Brent’s method is due to Richard Brent[1] and builds on an earlier algorithm by Theodorus Dekker.[2] Consequently, the method is also known as the Brent–Dekker method.

其内部组合二分, 割线, 逆抛物内插

精确线搜索

对精确的line search（线搜索），有一个重要的定理：

$$\nabla f(x_k+\alpha_kd_k)d_k=0$$

也就是下一个迭代点的梯度和当前搜索方向是垂直的.

证明如下

$f(\alpha_k)#是最小值, 那么有

$$f’(\alpha_k)=0=f’(x_k+\alpha_kd_k)d_k=\nabla f(x_k+\alpha_kd_k)^Td_k$$

非精确线搜索的两个准则

1. armijo-goldstein
2. wolfe-powell

mit armijo
对第一个准则

$$ f(x_k+\alpha_kd_k)\leq f(x_k)+\alpha_k \rho g_k^Td_k$$

$$f(x_k+\alpha_kd_k)\geq f(x_k)+\alpha_k (1-\rho) g_k^Td_k$$

首先这里注意一个充要条件

$$f(x_k+\alpha_kd_k)=f(x_k)+g_k^Td_k+o()$$

也就是说要下降, 必须$g_K^Td_k\leq0$

下面这张图很好的解释了armijo条件

可以看到armijo条件也有可能把极小值排除在外, 所以引入wolfe-powell条件

直观的wolfe条件就是下一个点处的梯度必须相比原来的点的斜率更平缓.

Trust Region

数值微分

1. 有限差分
2. 自动微分

trust region

Trust region methods are in some sense dual to line search methods: trust region methods first choose a step size (the size of the trust region) and then a step direction while line search methods first choose a step direction and then a step size.

信赖域方法相对线搜索方法的另一种方法, 可以说是一种对偶方法, 先确定步长(就是信頼域的规模), 然后确定方向, 而线搜索先确定方向, 在确定步长.LM方法就是信頼域方法.

In mathematics and computing, the Levenberg–Marquardt algorithm (LMA or just LM), also known as the damped least-squares (DLS) method, is used to solve non-linear least squares problems. These minimization problems arise especially in least squares curve fitting.

The LMA is used in many software applications for solving generic curve-fitting problems. However, as for many fitting algorithms, the LMA finds only a local minimum, which is not necessarily the global minimum. The LMA interpolates between the Gauss–Newton algorithm (GNA) and the method of gradient descent. The LMA is more robust than the GNA, which means that in many cases it finds a solution even if it starts very far off the final minimum. For well-behaved functions and reasonable starting parameters, the LMA tends to be a bit slower than the GNA. LMA can also be viewed as Gauss–Newton using a trust region approach.

用二次型来取代原来的函数, 但是求解位移的时候有一个信頼域限制, 这里半径利用预测下降量()和实际下降量的比值作为比值, 比值

注意NM算法是单纯形法!.

梯度法和牛顿法

梯度法只利用一阶信息, 利用平面逼近局部, 因为梯度是函数值变化最大的一个方向, 减去就可以不断下降.

牛顿法利用二阶信息, 利用曲面逼近局部. 因为牛顿法是在迭代点处用一个二次型函数代替原来的函数, 通过求这个二次型函数的极小值来作为下一次迭代点, 而二次新函数对x来说是一个曲面.
对目标函数二阶泰勒展开

$$f(x) = f(x^k)+(x-x^k)^Tg^k+\frac{1}{2}(x-x^k)^TF(x^k)(x-x^k)$$

其中令一阶导数为0,

$$x^{k+1}=x^k-F^{-1}(x^k)g^k$$

收敛性

梯度法如果采用最速下降法, 有全局收敛性, 对定步长, 需要步长小于一定值, 收敛阶数为1.

牛顿法的收敛性需要黑塞矩阵正定, 优点是当初始点离极值比较近的时候, 收敛性比较好, 对二次型函数, 一次迭代就可以收敛到极值, 收敛率阶数为2, 但是当离极值点比较远的时候, 并不一定收敛, 甚至反方向变化(黑塞矩阵奇异), 需要进行修正, 引入步长进行修正, 原理是牛顿法给出的改变方向是下降方向(因为黑塞矩阵大于0, 也就是正定矩阵, 方向一定是下降方向, 只要步长不太大.), 但是有可能步长导致不下降, 有三个缺点.

1. 要求黑塞矩阵
2. 要求逆
3. 初始点需要离极值点足够近(因为是用二次型函数近似原来的函数, 只有当足够近的时候, 这种近似才足够好)

LM修正-保证方向为下降方向

进一步由LM修正, 保证黑塞矩阵正定, 即保证每次得到的方向是下降方向, 给黑塞矩阵加单位矩阵就可以.

$$x^{k+1}=x^{k}-(F(x^k)+\mu_kI)^{-1} g^k$$

引入${\mu}_k$进行LM修正, 保证下降方向, 实际使用, 可以$0\rightarrow\infty$变化, 直到满足正定条件. 原理是加单位矩阵个数直接影响特征值, 进而可以使所有的特征值大于0, 也就变成了正定矩阵

最小二乘

最小二乘(线性或者非线性)可以直接用LM方法进行求解.

需要注意的是最小二乘拟合问题是$R^n\rightarrow R^n$问题, 求的不是梯度, 而是jacobian矩阵.

$$r=[r_1,…,r_m]^T$$

其中

$$r_i(x)=y_i-f(x_i)$$

即误差, 只需要$r$向量总的平方和最小即可.

那么其实雅克比矩阵每一列的函数是一样的, 因为只是数据不同.

目标函数的梯度用雅克比矩阵表示

$$\nabla f=2J(x)^Tr(x)$$

黑塞矩阵可以用雅克比矩阵进行表示

$$F(x)=2(J(x)^TJ(x)+S(x))$$

牛顿法表示为

$$x_{k+1}=x_k-(J(x)^TJ(x)+S(x))^{-1}J(x)^Tr(x)$$

其中$S(x)$包含$r$的二阶导数, 其中元素很小, 实际应用中直接忽略掉.就有了高斯-牛顿法

$$x_{k+1}=x_k-(J(x)^TJ(x))^{-1}J(x)^Tr(x)$$

进一步引入方向下降修正—-LM修正, 防止$J(x)^TJ(x)$不正定.

对高斯-牛顿法, 可以看到这个迭代公式中与线性最小二乘法解析公式基本是一致的, 不同的是把原来的数据$X,b$变成了$J(x),r(x)$, 分别是误差函数和其jacobian矩阵(本质是误差函数).

$$x_{k+1}=x_k-(J(x)^TJ(x)+\mu_kI)^{-1}J(x)^Tr(x)$$

共轭梯度

介于梯度法和牛顿法之间. 共轭的意思是关联, 基于两次搜索方向之间存在关联. 共轭向量之间是互不相关的向量. 相当于正交基.

计算上更好

1. 不用计算黑塞矩阵
2. 不用计算矩阵的逆

收敛性差点.

对于n维2次型问题, n次之内得到结果, 而牛顿法一步得到结果.

总的来说, 性能上, 牛顿法最好, 其次是共轭梯度, 最后是梯度法, 这里效率的关键是每次迭代的搜索方向.

原理是利用二次型的共轭方向作为搜索方向, 二次型的共轭方向定义为

$$d^TQd=0$$

共轭方向可以由很多, 并且线性无关, 个数$\leq n$

注意一个特性是每一步的搜索方向和下一步的梯度正交, 只要步长是最优步长(理论值需要计算梯度和黑塞矩阵).

由于每一步需要计算梯度和黑塞矩阵, 对于非线性函数, 比较难算, 把$\alpha$替换成一维搜索, 进一步利用修正方法计算黑塞矩阵, 只利用梯度信息.

总的来说, 利用迭代法求解极值得时候, 第一需要保证搜索方向是下降方向, 也就是黑塞矩阵为正定的, 第二保证具有下降特性, 就需要进行线性搜索.

最小二乘法

几何上, 理解为找到$b$在$A$上的正交投影, 因为当$m>n$时, $A$形成的子空间非常薄, 向量$b$在$A$的列向量形成的子空间的概率非常小, 而我们最小二乘要找的就是在$A$形成的空间中找到一个向量使得其到$b$向量的距离比其他向量都要小, 也就是正交投影向量, 注意当$m>n$时, $A$代表的空间为一个m维空间的一个n维子空间, 而当$m=n$的时候, $A$代表了整个n维空间, 所以$b$必然在其内部.

求解上, 把目标函数写成

$$f(x)=(Ax-b)^T(Ax-b)$$

这是一个二次型函数, 直接令其$\nabla f=0$就可以找到极小值对应的$x^ *$

进一步, 可以导出正交投影算子

$$P=I-A^T(A^TA)^{-1}A$$

进一步的, 对于$n>m$, 存在无数组解, 只能求最小范数解$||x||$, 解析解为

$$x^*=A^T(A^TA)^{-1}b$$

更一般的, 不需要分情况, 直接求出伪逆即可

$$x^*=A^+b$$

该结果有两个性质

1. 最小二乘对应的值最小.
2. $x^*$的范数最小.

loess

当优化函数是

$$L=(Xb-Y)^TD(Xb-Y)$$

其中$D$是$x_i,x_j$相似度的某种度量, 可以表示为内积, 和LR一样求导, 可以得到理论解

$$b=(X^TDX)^T(X^TDY)$$

这里$d(i,j)$可以看做一种核, 最常见的是高斯核

$$d(i,j)=e^{-\frac{||x_i-x_j||}{2\eta^2}}$$

可以看到是loess是LR的推广, 这里理论解都可以通过求导求解, 因为是凸问题.

伪牛顿法

主要包含

1. DFP
2. BFGS
3. l-bfgs

一个总结

伪牛顿法都是利用一个近似的正定矩阵去代替黑塞矩阵的逆.

其中BFGS是在DFP的基础上在增加的一些项.

每一迭代需要计算的量都是$s_k=x_{k+1}-x_k$,和$y_k=g_{k+1}-g_k$,和上一步的$H_k$

其中DFP算法如下

BFGS算法如下

l-bfgs算法如下

前两者需要保存n*n的矩阵,当这个矩阵维度比较大的时候占用内存比较大, 引入l-bfgs, limited-bfgs, 需要保存前$m$步的位置变化和梯度变化, 其他的都可以在此基础上计算, 算法如下

It is also possible to run BFGS using any of the L-BFGS algorithms by setting the parameter L to a very large number.

scipy的优化算法

l-bfgs一般使用的线性搜索算法是回溯线搜索.

In (unconstrained) minimization, a backtracking line search, a search scheme based on the Armijo–Goldstein condition, is a line search method to determine the maximum amount to move along a given search direction. It involves starting with a relatively large estimate of the step size for movement along the search direction, and iteratively shrinking the step size (i.e., “backtracking”) until a decrease of the objective function is observed that adequately corresponds to the decrease that is expected, based on the local gradient of the objective function.

算法步骤

全局算法

单纯形法

单纯性是指n维空间中由n+1个点组成的几何形状, 并且点满足

$$det[q_0,…,q_{n};ones(1,n+1)]\neq0$$

单纯性的迭代过程包括

1. 去掉最差的点后求重点.
2. 利用重心对最差的点进行反射.

这里反射根据反射点的优良性质又有很多不同的操作, 比如延伸(变成最好的点), 收缩(变成最坏的点), 压缩(收缩失败, 进行压缩).

对相等函数值的点, 采用平分决胜, 新的点的索引值更大(也就是越差).

模拟退火算法

一种随机搜索算法, 最主要的引入的概念是接受概率, 这个概率分布一般是基于gibbs分布, 这个分布的参数就是温度, 而温度是随迭代次数变化下降的, 所以称为退火.

$$p(k,f(z^k),f(x^k))=min(1,e^{-\frac{f(z^k)-f(x^k)}{T_k}})$$

同时这个概率分布的性质就是随着温度参数的下降, 接受坏点的概率也下降, 同时这个较差点离原来的点越远, 接受概率也越小.

常用语求解组合优化问题.(可行集有限, 但是

粒子群算法

也是随机搜索算法, 但是一次不是更新一个点, 而是一群点, 每一个点除了, 还有一个速度属性, 然后维护两个量, 一个是单个粒子历史最好位置向量, 一个是全局历史最好位置, 每次迭代都利用这两个量对每个例子的速度进行更新.

直观的, 维护四个量, 位置, 速度, 每个个体历史最好, 所有个体历史最好, 进行迭代时, 先利用两个历史最好量更新速度, 然后利用速度更新位置, 然后利用更新的位置更新两个历史量, 不断重复. 更新速度时候有一个惯性$w$, 一般取比1小点的数, $c_1,c_2$一般取2.$c_1$代表认知, $c_2$代表社会.

$$v_i^{k+1}=wx_i^k+c_1r_i^k{\circ(p_i^{k}-x_i^{k})+c_2s_i^k(g^k-x_i^k)}$$

$$x_i^{k+1}=x_i^k+v_i^{k+1}$$

遗传算法

也是基于种群的算法. p(0)到p(1)的过程, 每次迭代进行一次选择和进化

其中选择有两种模式, 一种是轮盘赌(归一化), 一种是竞标赛(随机选两个比较)

进化时开展交叉和变异, 首先确定父代, 然后进行交叉操作, 最简单的就是生成一个随机数形成交叉点, 直接把两个父代交叉点右边的全部交换, 也可以生成两个交叉建, 中间段互换. 接下来进行变异, 从配对池逐个抽个体进行变异, 一般变异的几率比较小, 所以影响没有交叉大.

这些全局算法的差异在

1. 不用计算导数
2. 每次迭代采用的都是随机操作
3. 利用一组点, 而不是个点

4. 编码(可采用实数编码)

   实数编码需要对应的交叉变异操作, 最简单的交叉可以直接取两个父代的平均值随机取代一个. 最简单的变异可以随机加一个均值为0的随机向量, 也可以利用约束集内任意一点的凸组合, 注意凸组合

   $$x’={\alpha}x+(1-\alpha)w$$

   其中$\alpha\in(0,1), w\in\Omega$, $\Omega$是约束集.

差分进化

线性规划

最简单的约束优化就是线性规划, 进一步最简单的凸规划, 就是凸二次规划, 两者都可以用单纯形法进行求解.

引入线性规划的对偶形式, 我们可以引出内点法(使用规模比较大的计算), 同样的对于凸二次规划, 我们可以由积极集法(positive-set)

有约束规划

引入KKT条件(和拉格朗日是一起的)

先讨论凸优化问题, (因为这类问题可以得到全局最优解)

什么是凸集?就是在凸组合这个操作下是封闭的.

凸函数?

$$\forall x,y\in\Omega, \alpha\in(0,1)$$

我们有

$$f(\alpha x+(1-\alpha)y)\leq \alpha f(x)+(1-\alpha)f(y)$$